#include <stdio.h>
#include <graphics.h>
#include <math.h>


double*solvematrix(double**M,double*P,int N)
{
    double*X=new double[N];
 
    int i, j, m; 
    double K, S = 0;
    double A[N][N];  
    double B[N]; 
    
    for (j = 0; j < N; j++) 
    {
      for (i = 0; i < N; i++) 
      {
        A[i][j] = *(((double*)M)+i*200+j);
      }
      B[j]=P[j];
    }

    for (i = 0; i <= N - 2; i++)
    {
      for (j = i + 1; j <= N - 1; j++)
      {
        K=A[i][i]/A[i][j];
        for (m = i; m <= N - 1; m++) 
        {
          A[m][j] = A[m][j] * K - A[m][i];
        }
        B[j] = B[j] * K - B[i];
      }
    }
    
    for (j = N - 1; j >= 0; j--) 
    {
      S = 0;
      for (i = N - 1; i > j; i--) 
      {
        S += A[i][j] * (*(X + i));
      }
      *(X + j) = (B[j] - S) / A[j][j];
    }
    
    return X;
}

int main() 
{
    int N = 4; 
            
    double A[200][200];
    double B[200];
    double*X=new double[N];

    double Norm;
    int Step = 0;
    double e = 0.01;
    int i, j = 0;
    
    initwindow(400, 300, "Norm(N)");
    setcolor(BLACK); 
    setbkcolor(WHITE);

    setfillstyle(1,RED);
    cleardevice();
    
    int x0 = 100, y0 = 200;   
    int x1 = 300, y1 = 100;   

    moveto(x0,y0);
    lineto(x0,y1);
    moveto(x0,y0);
    lineto(x1,y0);

    moveto(x0-3,y1+5);
    lineto(x0,y1);
    moveto(x0+3,y1+5);
    lineto(x0,y1);
    
    moveto(x1,y0);
    lineto(x1-5,y0-3);
    moveto(x1-5,y0+3);
    lineto(x1,y0);

    X = solvematrix((double**)A, (double*)B, N);
    
    for (i = 0; i <= N-1; i++)
    {
      printf("x%d = %f; ",i + 1,*(X+i));
    }
    printf("\n");
    
    N = 1;
    while(true==true)
    {
      N++;
      double * X1 = new double[N];
      printf("N = %d\n\n", N);
      for (j = 0; j <= N-1; j++)
      {
        for (i = 0; i <= N-1; i++)
        {
          if(i==j)  
          {
            A[i][j] = 1 + e;
          }
          else if(i<j)
          {
            A[i][j]=e;
          }
          else if(i > j)
          {
            A[i][j] = 1 - e;
          }
        }
        if(j==N-1)
        {
          B[j] = 1 + e;
        }
        else
        {
          B[j] = -1;
        }
      }   
      
      for (j = 0; j <= N-1; j++) 
      {
        for (i = 0; i <= N-1; i++) 
        {
          printf("%f ", A[i][j]);
        }
          printf("-- %f\n", B[j]);
      }   
      printf("\n"    );

      X1 = solvematrix((double**)A, (double*)B,N);
      for (i = 0; i <= N-1; i++)
      {
        printf("x%d = %.15f\n",i + 1,*(X1+i));
        Norm+=X1[i]*X1[i];
      }
      printf("\n");
        
      Norm = pow(Norm, 0.5) / N;
      printf("Norm = %.15f\n\n\n", Norm);
      fillellipse(101 + Step*5, 200 - (int)(Norm*100), 1, 1);
      Step++;
      delete X1;
      if(Step>=16)
      {
        break;   
      }
    }  
    
    getchar();
    closegraph();
	return 0;
}
